Parameter Identification for a One-dimensional Blood Flow Model

The purpose of this work is to use a variational method to identify some of the parameters of one-dimensional models for blood flow in arteries. These parameters can be fit to approach as much as possible some data coming from experimental measurements or from numerical simulations performed using more complex models. A nonlinear least squares approach to parameter estimation was taken, based on the optimization of a cost function. The resolution of such an optimization problem generally requires the efficient and accurate computation of the gradient of the cost function with respect to the parameters. This gradient is computed analytically when the one-dimensional hyperbolic model is discretized with a second order Taylor-Galerkin scheme. An adjoint approach was used. Some preliminary numerical tests are shown. In these simulations, we mainly focused on determining a parameter that is linked to the mechanical properties of the arterial walls, the compliance. The synthetic data we used to estimate the parameter were obtained from a numerical computation performed with a more accurate model: a three-dimensional fluid-structure interaction model. The first results seem to be promising. In particular, it is worth noticing that the estimated compliance which gives the best fit is quite different from the values that are commonly used in practice. Résumé. Le but de ce travail est d’identifier certains des paramètres existant dans des modèles 1-d d’écoulement sanguin dans des artères. Ces paramètres peuvent permettre d’approcher autant que possible des configurations géométriques réalistes ou des données expérimentales. Une approche de l’estimation de paramètres par moindres carrés non-linéaires a été adoptée, basée sur l’optimisation d’une certaine fonction coût. La résolution d’un tel problème de minimisation requiert le calcul efficace et précis du gradient de la fonction coût par rapport aux paramètres. Le gradient est discrétisé analytiquement dans le cas d’une discrétisation du modèle hyperbolique 1-d par le schema de TaylorGalerkin. Une approche par l’état adjoint a été employée. Des premiers résultats numériques sont fournis. Pour ces simulations, nous nous sommes concentrés sur la détermination d’un paramètre lié aux propriétés mécaniques de la paroi artérielle. Les données synthétiques utilisées pour l’estimation de ce paramètre ont été obtenues à partir d’un modèle beaucoup plus raffiné : un modèle 3-d d’interaction fluide-structure. Les résultats semblent intéressants car le paramètre estimé est assez différent de ce à quoi on s’attendrait a priori. ∗ This study was supported by RTN-Project “HaeMOdel”, project no HPRN-CT-2002-00270. 1 MOX, Dpto di Matematica, Politecnico di Milano, Via Bonardi 9, 20133 Milano, Italy ; e-mail: [email protected] & [email protected] 2 Estime, Inria Rocquencourt, BP 105, 78153 Le Chesnay, France ; e-mail: [email protected] 3 Bang, Inria Rocquencourt, BP 105, 78153 Le Chesnay, France ; e-mail: [email protected] 4 Reo, Inria Rocquencourt, BP 105, 78153 Le Chesnay, France ; e-mail: [email protected] c © EDP Sciences, SMAI 2005 Article published by EDP Sciences and available at http://www.edpsciences.org/proc or http://dx.doi.org/10.1051/proc:2005014 PARAMETER IDENTIFICATION FOR A ONE-DIMENSIONAL BLOOD FLOW MODEL 175 Introduction We focus in this study on the parameter estimation of a 1-d blood flow model, [11, 19]: ∂A ∂t + ∂Q ∂z = 0, ∂Q ∂t + ∂ ∂z ( αQ A ) + A ρ ∂P ∂z +Kr ( Q A ) = 0, where the pressure P , the area A and the flux Q are the unknowns of the problem. We denoted by z the abscissa, t the time, ρ the density of the blood. Two parameters, the Coriolis coefficient α and the friction parameter Kr, are introduced in this model. This system is closed with a wall displacement law of the form P (t, z)− Pext = β̃ ( A −A 0 ) , that links the pressure and the area. We introduced the external pressure Pext, the area at rest A0 and a coefficient β̃ that takes into account the mechanical characteristics of the arterial wall. The parameters (α, β̃,Kr, A0) used in this model are related to physiological data or to the velocity profile. Thus our aim is to identify some of these parameters. Motivations. Two main objectives can be thought of to motivate the parameter estimation in 1-d models. The first goal may have interesting clinical applications. Knowing some non-invasive clinical data measured on a patient, one would like to retrieve the actual physiological or mechanical constants of this patient. For instance, it is possible to measure unintrusively the mean fluxes and areas as a function of time at two or three different sections of an artery. From these data, one would like to identify the mechanical properties of the arterial wall: we could thus obtain the compliance of the wall and the pressure in a totally non-invasive manner. The second objective is consists in making a coarsening of models. One is now able to solve the full 3-d fluid structure interaction problem on real geometries, coming from real patients. However the resolution of this problem is quite expensive and this complex model cannot probably be afforded for an intensive numerical study requiring lots of resolutions. This can occur when one wishes to modify the configuration of the flux or the boundary conditions for instance. In this case, a single resolution of the accurate but expensive 3-d model could provide data, such as the flux and the area for all sections of the mesh. Then one can estimate the 1-d parameters from the 3-d data; this would allow to make use of the cheap 1-d model for intensive computations, but using the parameters that take into account data coming from a real geometry and a physically more detailed model. This multiscale approach has already been used in a medical application: some 3-d NavierStokes (without compliant walls) were performed to obtain a simple numerical/experimental law that was used in ulterior 0-d simulations, [16, 17]. The difference is that in this study, we base our parameter identification on sound mathematical tools, and we focus on 1-d models that provide a more accurate description of wave propagation in the large arteries. One of the conclusions of this work is that the coefficient estimated by solving the inverse problem is quite different from the coefficient one would have chosen a priori (an a priori expression for β̃ is provided in Section 1). Thus this approach could provide 1-d models that are suitable for a 3-d–1-d coupling in multiscale computations, see [8]. This illustrates the relevance of our approach. To conclude on the motivations, one can either aim at estimating physical parameters from experimental unintrusive measurements, or at having a cheap 1-d model to be as close as possible to an expensive 3-d model, in order to make realistic configuration studies using this cheap model. We note here that the methodology remains identical for both objectives, the only difference being the expression of the measurement operator. Methodology. We followed a standard nonlinear least squares approach to parameter estimation, [3–5]. It is based on the optimization of an appropriate cost function. The resolution of such a minimization problem generally requires the efficient and accurate computation of the gradient of the cost function with respect to the parameters, [4]. We discretized the 1-d model with a second order Taylor-Galerkin scheme, [1]. We computed analytically the gradient of the discrete cost function, using an adjoint approach. The adjoint problem obtained 176 VINCENT MARTIN, FRANCOIS CLÉMENT, ASTRID DECOENE AND JEAN-FRÉDÉRIC GERBEAU is, as expected, a linear 1-d hyperbolic system, but has nonstandard discretization and boundary conditions, that are due to the differentiation of the Taylor-Galerkin scheme. A previous attempt was made in [14] to estimate the elasticity of the arteries. However, instead of differentiating the discrete equations as we do here, the author differentiated first the continuous equations and then discretized the continuous adjoint problem thus obtained. We believe this is a possible reason for the very slow convergence he reached in the minimization process (about 1500 iterations for three parameters). Finally, this approach allows to make a sensitivity analysis, [13, 15], that can provide information on the relevant parameters to estimate, and on the type of measurements to perform. But this goes beyond the scope of this study. Results. The analytical discrete gradient was implemented and validated by comparison with finite differences approximations. The adjoint problem was not computed by automatic code differentiation, but directly implemented. Thus we could easily control the memory required by the gradient code. We used a constrained optimization code based on a quasi-Newton method with active constraints. We present some preliminary numerical results. In these numerical simulations, we mainly focused on determining the parameter β̃ that is linked to the mechanical properties, i.e. the compliance, of the arterial walls. The synthetic data we used to estimate the parameter were obtained from a numerical computation performed with a 3-d fluid structure interaction model. We first used as data the values of the areas and fluxes at only two or three points of the domain (boundaries plus maybe the middle point). Although the data at two points do not seem to be enough to find a stable value, it seems that with three points, one can obtain a β̃ relatively stable (i.e. it is little changed when the estimation is made with all available spatial data). In the second numerical tests, we used all data available from the 3-d computation. These first numerical results seem promising and should be followed by further developments. In Section 1, we present briefly the continuous 1-d model, that is derived in Section 2 with the Taylor-Galerkin scheme. In Section 3, the gradient of the least squares cost function is computed with the adjoint approach. Numerical results are presented in Section 4, and some conclusions and perspectives in Section 5. 1. Direct model: 1-d blood flow model We present in this section a 1-d blood flow model based on the works in [10, 20]. See also [19]. It is a 1-d vectorial hyperbolic problem, with a 2 × 2 flux matrix that admits two real eigenvalues with opposite signs under physiological conditions. We leave the problem of the parameterization and of the measurements for the next sections, see Section 3.1. 1.1. Continuous blood flow model Let Ω = (0, L) be a 1-d domain of length L > 0. Let I = (0, Tf), with Tf > 0, the time interval of simulation. The continuous system of equations reads, for the abscissa z ∈ Ω and the time t ∈ I ∂A ∂t + ∂Q ∂z = 0, z ∈ Ω, t ∈ I, ∂Q ∂t + ∂ ∂z ( αQ A ) + A ρ ∂P ∂z +Kr ( Q A ) = 0, z ∈ Ω, t ∈ I, (1) where U = [A,Q] is the vectorial unknown of the problem, made of the area A and of the flux Q. The pressure P is an intermediary unknown. We denoted by ρ the blood density that we assume perfeclty known in this study. Two parameters, the Coriolis coefficient α and the friction parameter Kr, are introduced when deriving this model. This system of partial differential equations is completed with an initial condition U(0, z) = U0(z), z ∈ Ω, (2) PARAMETER IDENTIFICATION FOR A ONE-DIMENSIONAL BLOOD FLOW MODEL 177 and some adapted boundary conditions φ0(U(t, 0),p) = q0(t), t ∈ I, φL(U(t, L),p) = qL(t), t ∈ I . (3) In equations (2) and (3), U0, q0 and qL are given initial and boundary data. The definition of the real-valued boundary functions φ0, φL is discussed in Section 2.2. The system (1) is closed with a wall displacement law of the form, see [21] for instance and the reference therein, P (t, z)− Pext = ψ(A;A0,β) = β0 [( A A0 )β1 − 1 ] , (4) where β = (β0, β1) is a pair of positive real parameters. The power coefficient β1 is often taken equal to 1/2, which means that the pressure difference is proportional to the wall displacement, and, in this case, a linear elastic law can provide an expression for β0, if the mechanical properties of the arterial wall are known a priori, β1 = 1/2, β0 = √ πh0,wE √ A0(1− ν2) , (5) where h0,w is the wall thickness, E is the wall Young modulus, and ν is the Poisson coefficient. We can reformulate the wall displacement law in the following way: P (t, z)− Pext = ψ(A;A0, β) = 2ρβ [ A −A 0 ] , with β = β0 2ρA0 = √ πh0,wE 2ρA0(1− ν2) . (6) We introduce the following quantity c = c(A;A0,β) = √ A ρ ∂ψ ∂A = √ β0β1 ρ ( A A0 )β1 , (7) which has the dimension of a velocity and is related to the speed of propagation of simple waves along the tube. We also introduce the integral of the square of the celerity c with respect to the area